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Abstract 

Molecular dynamics simulations are performed to investigate heterogeneous dynamics in amor- 
phous glassy materials under oscillatory shear strain. We consider three-dimensional binary 
Lennard-Jones mixture well below the glass transition temperature. The structural relaxation 
and dynamical heterogeneity are quantified by means of the self-overlap order parameter and the 
dynamic susceptibility. We found that at sufficiently small strain amplitudes, the mean square 
displacement exhibits a broad sub-diffusive plateau and the system undergoes nearly reversible 
deformation over about 10^ cycles. Upon increasing strain amplitude, the transition to the diffu- 
sive regime occurs at shorter time intervals and the relaxation process involves intermittent bursts 
of large particle displacements. The detailed analysis of particle hopping dynamics and the dy- 
namic susceptibility indicates that mobile particles aggregate into clusters whose sizes increase at 
larger strain amplitudes. Finally, the correlation between particle mobilities in consecutive time 
intervals demonstrates that dynamic facilitation becomes increasingly pronounced at larger strain 
amplitudes. 

PACS numbers: 61.43.-j, 66.30.Pa, 62.20.F-, 83.10.Rs 
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I. INTRODUCTION 



Understanding the relationship between atomic structure and mechanical properties in 
amorphous materials is important in many current applications and emerging technolo- 
gies ij. In contrast to crystalline solids, where the plastic deformation is governed by dislo- 
cations, it was originally found that the plastic activity in amorphous materials is controlled 
by the localized shear transformations [2], which were more recently studied using com- 
puter simulations 3-6| and directly visualized in experiments on colloidal glasses ^jlsl and 
foams [9^. Under the applied strain, the deformation of amorphous materials is determined 
by the cooperative organization of irreversible rearrangements of small clusters of particles, 



which could be triggered by the nonlocal redistribution of elastic stress 



10|-ll3|- The se- 



quence of such plastic events can lead to an avalanche process characterized by a power-law 
scaling of the average stress or energy drops with the system size [l4, [isl. Another notable 
examples of systems with intermittent, spatiotemporal heterogeneous dynamics include the 
Barkhausen crackling noise in magnets subject to an external magnetic field 16|] and the 



17, 



18|- 



driven block-spring model in the theory of self-organized criticality 

In recent years, computer simulations have become an increasingly important tool for 
studying slow particle dynamics in molecular liquids near the glass transition at thermal 
equilibrium [19]. A common observation is that the mean square displacement of individual 
particles is reduced upon approaching the glass transition temperature from above, which 
also results in broadening of the sub-diffusive plateau that separates the ballistic and diffu- 
sive regimes. More importantly, however, is that the particle mobility in the sub-diffusive 
regime ea„ be signifieantly differe,. from the average value; a„d. in addition, the partides 
with similar mobility become spatially correlated, thus leading to dynamic heterogeneity [2C 



22|. The spatial fluctuations of mobile regions are characterized by the four-point dynamic 



correlation function, and they can be efficiently measured by computing the variance of 



the self-overlap order parameter, or the dynamic susceptibility [23|, |24] . While no obvious 
changes in spatial correlations of particle positions are detected near the glass transition, 
the dynamic susceptibility provides an estimate of the number of particles involved in the 
correlated motion. It is now well recognized that when a liquid is cooled toward the glass 
transition temperature, the peak value of the dynamic susceptibility increases, indicating 



that dynamics becomes spatially increasingly correlated [19|. More recently, it was shown 
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that in the presence of steady shear flow, the dynamics of supercooled hquids is more ho- 
mogeneous as the shear flow reduces the dynamic correlation length and the lifetime of 
dynamical heterogeneity 25 1. 

The microscopic mechanism of structural relaxation in glassy materials is governed by 
spatially extended domains of fast moving particles that can be identifled from the analy- 
sis of individual particle trajectories. Hopping particle dynamics was recently studied in a 



number of systems, i.e., binary mixtures below the glass transition 
at thermal equilibrium dense granular media 



28 



26|, supercooled liquids 



29| . and actively deformed polymer 



glasses jsot. In each particle trajectory was decomposed into a series of segments, 

where motion takes place inside a cage, separated by fast cage jumps. Of particular impor- 
tance is the analysis of cage jumps and their spatio-temporal clusterization performed in the 
cyclic shear 



and the fluidized bed 



29 



31| experiments of two-dimensional granular me- 



dia. In both experiments, it was found that the major structural relaxation events are well 



correlated with the bursts of cage jumps. Furthermore, these cage _ 



into clusters, whose sizes are approximately power-law distributed 28 



umps tend to aggregate 



291]. In turn, several 



clusters might dynamically facilitate each other and form avalanches, which propagate along 



the soft regions of the system 



27H29|. One of the motivations of the present study is to 



examine the collective motion of particles and dynamic facilitation in periodically deformed 
three-dimensional amorphous materials. 

In this paper, molecular dynamics simulations are employed to study the relaxation dy- 
namics in the Kob-Andersen Lennard- Jones binary mixture model at a flnite temperature 
well below the glass transition. The three-dimensional system is periodically strained over 
many cycles, probing regions with low instability thresholds, which leads to intermittent 
localized rearrangements of particles. We flnd that at sufficiently large strain amplitudes 
the particle dynamics is purely diffusive, while at lower amplitudes, the mean square dis- 
placement develops an extended sub-diffusive plateau. The analysis of the dynamical sus- 
ceptibility and particle hopping dynamics reveals the spatial heterogeneity of structural 
relaxation. 

The rest of the paper is organized as follows. The details of molecular dynamics sim- 
ulations are described in the next section. The results for the particle diffusion, hopping 
dynamics, and microstructure of clusters of mobile particles, as well as the analysis of the 
two and four-point correlation functions are presented in Sec. lIIII The conclusions are given 
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in the last section. 



II. MOLECULAR DYNAMICS SIMULATION MODEL 

The three-dimensional amorphous material is modeled as the Kob-Andersen binary 
(80:20) Lennard- Jones mixture with non-additive interaction parameters that prevent crys- 



tallization 32|]. The snapshot of the equilibrated system which consists of Np = 2940 
particles is presented in Fig.[T] In this model, any two particles a, P = A, B interact via the 
pairwise Lennard- Jones (LJ) potential 



(1) 



where Eaa = 1-0, Eab = 1-5, Ebb = 0.5, aAs = 0.8, cbb = 0.88, and = m^. The cutoff 



radius is taken to be twice the minimum position of the LJ potential r^Q/? = 2.245 cTq,/? 



33|, 



34l |. In what follows, the units of length, mass, and energy are set to be a = aAA, = niA, 
and e = Eaa, and, correspondingly, the unit of time is defined as r = o^fmfe. The 
equations of motion were solved numerically using the fifth-order Gear predictor-corrector 



algorithm |35( with a time step tsiuD = 0.005 r 

All simulations were performed at a constant volume with the total density p = Pa + Ps = 
1.2 and temperature O.le/fcs, where ks is the Boltzmann constant. This temperature 
is well below the value 0.45 e/fc^ at which the computer glass transition is detected 32 |. 
The constant temperature was maintained by rescaling the velocity component in the y di- 
rection (perpendicular to the plane of shear). As indicated in Fig.[Tl the system dimensions 
are measured L^. = 12.81 a, Lj^ = 14.79cr, and Lz = l2.9Aa. In order to simulate homoge 



neous, time-dependent shear strain, the Lees-Edwards periodic boundary conditions [35 1 
were implemented with the SLLOD equations of motion [36]. It should be mentioned that 
in contrast to the boundary-driven shear algorithms, the spatially homogeneous shear strain 
prevents the formation of shear bands [s^]. 

The time-periodic shear strain was imposed (in the xz plane) by varying the strain rate 
as a function of time 7(t) = 7ocos(wt), where u is the oscillation frequency and 70 is the 
strain rate amplitude. We define the strain amplitude as a ratio of the strain rate amplitude 
over the frequency, i.e., 70 = ■jo/oj- For the results reported in this paper, the oscillation 
frequency was set ut = 0.02, with the corresponding period T = 27t/u = 314.16 r; and the 
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strain amplitude was varied in the range 70 ^ 0.08. 

The system was first equihbrated for about 5 x 10^ MD steps at a constant volume and 
temperature 1.2e/kB in the absence of shear, and then gradually quenched to the final 
temperature O.le/fc^ with steps of 0.1 e/A;^. After the oscillatory shear strain was applied, 
the first 2 x lO'^ MD steps were discarded to avoid quench-rate and aging effects. During the 
oscillatory motion, the measurements of particle positions were taken every back and forth 
cycle when strain is zero. The data were accumulated over 12, 000 cycles (about 7.5 x 10^ 
MD steps) at each strain amplitude, and the post-processing analysis of particle trajectories 
was performed in six independent systems. 



III. RESULTS 



At the studied temperature and density, an equilibrated model glass in the absence of 
deformation is characterized by the amorphous liquid-like molecular structure where most 
of the atoms remain in cages formed by their neighbors on the time scale accessible to 
computer simulations 32|. A typical steady shear stress-strain response involves an elastic 
deformation at strains below a few percent and a yield stress that depends on the physical 



aging and strain rate [3J, |38|. During the elastic and plastic deformations, the atoms can 
undergo non-affine, irreversible displacements, which, depending on the strain, might form 
cascades spanning a considerable fraction of the system jsol. Thus, instead of tempera- 
ture, the consecutive irreversible displacements of atoms are governed by the strain rate 
as a control parameter 



40]. In the present study, the amorphous material is periodically 
deformed, and the particle positions are saved every cycle at zero strain. In such a setup, 
the affine deformation field, which is present at steady shear strain, is zero; and, therefore, 
any irreversible particle rearrangements will contribute to the structural relaxation of the 
material. 

The mean square displacement (MSD) averaged over both A and B particles is plotted in 
Fig.[2]as a function of time for the frequency ut = 0.02 and various strain amplitudes. To 
compute the displacement, the position of the system center of mass was subtracted from 
the position of each particle. It can be observed in Fig.[2]that at small strain amplitudes, 
7o ^ 0.06, the MSD curves exhibit a broad sub-diffusive plateau, which becomes more 
pronounced at smaller strain amplitudes, and a gradual crossover to diffusive motion. In 
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contrast, at larger strain amplitudes, 70 ^ 0.07, the sub-diffusive regime is absent, and the 
slope of the MSD curves becomes equal to one at times t > 10 T as indicated by the straight 
dashed line in Fig.|2J At the largest time interval t = 1.2 x lO^'T, the particle displacement 
is still about the cage size for the strain amplitude 70 = 0.02, which implies that during the 
periodic deformation the system dynamics is nearly reversible. Note also that the ballistic 
regime is not observed in any of the MSD curves as it occurs at times much smaller than 
the oscillation period T = 314.16 r. 

The appearance of the extended sub-diffusive plateau in the MSD curves reported in 
Fig. [2] suggests that the particle dynamics might be spatially heterogeneous on the length 
scales of about the cage size. The structural relaxation in amorphous materials is commonly 
quantified via the self-correlation function, which is defined as follows: 

Q,(M)^i.E-p(-^). (2) 

i=l 

where Arj(t) = rj(to + t) — rj(to) is the displacement vector of the zth particle, t is the lag 
time, and a is the probed length scale j41|. In turn, the extent of dynamical heterogeneity 
is measured by the four-point correlation function, or the dynamical susceptibility, which is 
computed as the variance of the self-correlation function: 

X4(a, t) = Np [(Q,(a, tf) - {Q,{a, t)f] , (3) 



where the brackets (■) denote averaging over all particles and initial times 23j. At a given 



time lag, the correlation function X4(a,t) provides an estimate of the number of particles 



involved in a cooperative displacement over the length scale a 



19l |. At some intermediate 



time and length scales, the function ^4(0, t) usually displays a maximum indicating the 



23| 



largest spatial correlation between localized particles 

The time dependence of the self-correlation function Qs{a,t) is illustrated in Fig.[3]when 
LOT = 0.02 and the parameter a is slightly larger than the cage size, i.e., a = 0. 12 cr. As is 
evident, the correlation function Qs{ci, t) decays faster at larger strain amplitudes. Note that 
at smaller strain amplitudes 70 ^ 0.05, the system is not fully relaxed even at the largest 
time interval t = 1.2 x 10"^ T. On the other hand, at the smallest time interval t = T, the 
function Qs{a, t) is less than 1.0 because of the thermal vibrations inside a cage. Further, the 
dynamic susceptibility x^i^^y t) is shown in Fig.lHfor the same parameters ut = 0.02 and a = 
0.12 a. It can be observed that the correlation function X4('^) t) exhibits a pronounced peak, 



whose magnitude increases at larger strain amplitudes, indicating progressively larger size 
of dynamically correlated regions. Upon increasing strain amplitude, the peak is displaced 
to smaller times, which is consistent with the onset of the diffusive regime in MSD curves 
reported in Fig.|2J Finally, assuming that the correlated regions are compact, the dynamic 
correlation length ^4 can be estimated from the peak value of X4(a,t) at a = 0.12 a. The 
inset in Fig.H] shows ^4 as a function of the stain amplitude. It is seen that the data are 
well fitted by the power-law function with the exponent of about 0.9 (straight red line in 
Fig.H]). One should keep in mind, however, that the dynamic susceptibility is an averaged 
quantity, which measures mean square fluctuations of the number of mobile particles, and, 
thus, it does not completely describe the microscopic mechanism of structural relaxation. 
We next perform a more detailed analysis of the particle hopping dynamics and the local 
microstructure of clusters of mobile particles. 

At sufficiently low temperature and small strain, a typical particle trajectory in a glassy 
material consists of rapid hopping events separated by the rattling motion within a cage. 
Hence, the hopping dynamics is controlled by the cage-to-cage jumps, which, in practice, can 
be identified by a numerical algorithm recently introduced by Candelier et al. [28]. In essence, 
this algorithm is based on the spatial separation of two consecutive segments of a particle 
trajectory. More specifically, the measure of the distance separating two segments is defined 
by the product of the root mean square distances between all points within the segments to 
the center of mass of the other segment 281]. Furthermore, the effective distance between 
two segments is normalized by a factor that counterbalances large fiuctuations arising from 
averaging over short segments. During the iterative procedure, the cage jumps are detected 
if the effective distance is greater than the typical cage size; and the whole trajectory is 
consecutively divided into a number of segments where the particle motion takes place inside 
a cage. This algorithm was successively applied to identify cage jumps in two-dimensional 



granular systems under cyclic loading |28|], in the fluidized bed experiment [29|, and in 
supercooled liquids at mechanical equilibrium {2^]. 

In the present study, the cage detection algorithm was used to analyze three-dimensional 
trajectories of individual particles as follows. First, we choose a subset of points in the par- 
ticle trajectory, divide the subset in two adjacent segments, and then compute the effective 
distance separating these two segments. If the effective distance between any two adjacent 
segments within the subset is less than the cage size, r^. = 0.1 a, then the particle was con- 
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sidered being in the cage during the time interval defined by the subset. For every particle, 
the procedure was repeated for all time intervals greater than 10 T and less than 100 T and 
all initial times. As a result, all particle trajectories were decomposed into successive cages 
separated by cage jumps, which typically consist of several consecutive points each. A vi- 
sual examination of the trajectories revealed two types of jumps; namely, reversible, where 
a particle jumps back and forth between the averaged positions, and irreversible, where a 



particle permanently escapes its cage, 
studies of glassy systems at equilibrium 



his is consistent with the results of previous MD 



26 



42 



43|. 



FigureE] shows the total number of particles undergoing cage jumps as a function of 
time for three representative cases at strain amplitudes 70 = 0.02, 0.04, and 0.06. It is 
clearly observed that the periodic deformation generates a heterogeneous temporal response 
characterized by intermittent bursts of large particle displacements. It is apparent that the 
amplitude of the bursts and frequency of their occurrence increase at larger strain ampli- 
tudes. During the time intervals between the bursts, we also detect a finite number of cage 
jumps that are assisted by thermal activation. Following tradition, the frequency spectrum 
of the data series in Fig.[5]was determined by computing their discrete Fourier transform. 
Within the reported time interval, the power spectrum at each strain amplitude exhibits 
a power-law decay with the exponent of about two (not shown), which is indicative of a 
simple Brownian noise. This is in contrast with the inverse frequency spectrum of the so- 
called flicker noise found in many complex systems that are characterized by scale-invariant 



ava. 
ity 



anche-like processes and described by the phenomenon known as self-organized critical- 



17 



44|. 



A more direct evidence of spatial heterogeneity can be obtained from visualization of 
instantaneous positions of mobile particles. Snapshots of mobile particle positions during 
intermittent bursts are presented in Fig. [6] for different strain amplitudes. It is clearly seen 
that the particles undergoing cage jumps mostly aggregate into clusters whose sizes increase 
at larger strain amplitudes. A number of previous studies of heterogeneous dynamics in 
glassy materials have demonstrated that the sizes of clusters of mobile particles are power- 



law distributed 



28 



29 



44j |. As shown in Fig.[5l the structural relaxation process involves only 
a few large-scale cooperative clusters during the time interval 10^ T, and, as a consequence, 
we find that the probability distribution of cluster sizes of more than about 20 particles is 
subject to large statistical uncertainty (not shown). 
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It was recently suggested that dynamic facilitation might be one of the important mech- 
anisms leading to spatio-temporal heterogeneity in glassy materials 1^5^. In the context of 
kinetically constrained models, dynamic facilitation can be quantified either by the mobility 
transfer function between highly mobile regions and nearby regions that were previously 
mobile, or by the facilitation volume, which measures the spatial extent of mobile regions 
initiated by localized excitations 46|. In MD simulations of glass-forming liquids, the mo- 
bility transfer function was computed for mobile particles near their neighbors that were 



previously mobile [471, |48| . In particular, it was demonstrated that mobility propagates con- 



tinuously through the system, and dynamic facilitation becomes increasingly pronounced 



upon supercooling 



47 



48|. 



In the present study, two measures of dynamic facilitation were considered based on the 
results of the cage detection algorithm applied to individual particle trajectories. First, 
it was determined whether a particle was mobile at a given time step and it remained 
immobile during the preceding time interval At. Next, we checked if the particle had at 
least one mobile nearest neighbor during the time interval A t. Figure[7] shows the ratio 
of dynamically facilitated mobile particles and the total number of particles that become 
mobile after At/T cycles. It can be seen that the ratio Nf/Nfot increases with increasing 
strain amplitude, implying that dynamic facilitation plays a more important role at larger 
strain amplitudes. Note also that at the strain amplitude 70 = 0.06 and At/T > 10^, nearly 
all particles undergo cage jumps after being in contact with mobile regions. In contrast, at 
the strain amplitude 70 = 0.02, the ratio Nf/Ntot appears to saturate at about 0.65, which 
characterizes the relaxation dynamics that involves single particles undergoing reversible 
jumps and rearrangement of small clusters of particles [e.g., see Fig.pfa)]. 



Similar to the analysis presented in the previous MD studies |47|, |48| , we also computed 
the correlation between particle mobilities in back-to-back time intervals of equal duration 
A t. Namely, we only selected particles that were mobile at least once during the time inter- 
val A t but always immobile during the preceding time interval A t. Then, the dynamically 
facilitated mobile particles were identified if there was at least one mobile nearest neighbor 
during the preceding time interval At. The results are shown in the inset of Fig. [71 Al- 
though the data are somewhat noisy, the trend is clear; an increasingly larger fraction of 
mobile particles are dynamically facilitated at larger strain amplitudes. Thus, regardless of 
the definition, the results of numerical simulations indicate that, as the strain amplitude 
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increases, there is a higher probabihty to find a mobile particle that was previously located 
near mobile regions. 

IV. CONCLUSIONS 

In this paper, we performed molecular dynamics simulations to study heterogeneous re- 
laxation dynamics in an amorphous material under time-periodic shear strain deformation. 
The three-dimensional amorphous material was modeled as the binary Lennard- Jones mix- 
ture at a temperature well below the glass transition. During the oscillatory deformation, 
the particle positions were stored every cycle when the net strain is zero. We found that at 
small strain amplitudes, the mean square displacement develops an extended sub-diffusive 
plateau followed by the diffusive regime; whereas at larger amplitudes only the diffusive 
regime is present at the reported time scales. 

The structural relaxation was described by the decay of the self-overlap correlation func- 
tion, which indicated that at small strain amplitudes the system dynamics is nearly reversible 
over about 10^ cycles, while at strain amplitudes above a few percent, almost all particles un- 
dergo irreversible displacements and escape their cages. With increasing strain amplitude, 
the dynamic susceptibility exhibits a pronounced peak at intermediate time and length 
scales, and the magnitude of the peak increases at larger strain amplitudes, indicating pro- 
gressively larger size of dynamically correlated regions. Furthermore, the detailed analysis 
of particle hopping dynamics revealed that the periodic deformation generates a heteroge- 
neous temporal response characterized by intermittent bursts of large particle displacements. 
Lastly, our numerical simulations have shown that dynamic facilitation of mobile particles 
becomes increasingly important as the strain amplitude increases. 
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FIG. 1: (Color online) A snapshot of instantaneous positions of particles A (blue circles) and B 
(red circles) at mechanical equilibrium. The particle sizes are not drawn to scale. During the 
oscillatory motion, the periodic shear strain was applied in the xz plane (indicated by the arrows). 
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FIG. 2: (Color online) The mean square displacement of A and B particles as a function of time 
for the oscillation frequency (jjt = 0.02 and period T = 2tt/oj = 314.16 r. The strain amplitudes 
from bottom to top are 70 = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, and 0.08. The dashed line with 
unit slope is plotted for reference. 
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FIG. 3: (Color online) The self-correlation function Qs{a,t) defined by Eq. ([2]) for the oscillation 
frequency ujt = 0.02 and period T = 211/00 = 314.16 r. The probed length scale is a = 0.12 o". The 
strain amplitudes from top to bottom are 70 = 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, and 0.08. 
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FIG. 4: (Color online) The dynamic susceptibility X4('^)^) foi^ the oscillation frequency ojt = 0.02 
and a = 0.12 a. The strain amplitudes from bottom to top are 70 = 0.02, 0.03, 0.04, 0.05, 0.06. 
The oscillation period is T = 314.16 r. The inset shows the dynamic correlation length ^4 = 
[X4"^^(i)]^^'^ as a function of the strain amplitude 70. The red line with a slope 0.9 is the best fit 
to the data. 
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FIG. 5: (Color online) The number of particles undergoing cage jumps as a function of time for the 
oscillation frequency ujt = 0.02, period T = 27r/a; = 314.16 r, and strain amplitudes (a) 70 = 0.02, 
(b) 7o = 0.04, and (c) 70 = 0.06. Note that the vertical scale is different in the panel (c). 
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FIG. 6: (Color online) Typical clusters of mobile particles A (blue circles) and B (red circles) 
for the oscillation frequency cot = 0.02 and strain amplitudes (a) 70 = 0.02, (b) 70 = 0.03, (c) 
70 = 0.04, and (d) 70 = 0.05. 
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FIG. 7: (Color online) The ratio of dynamically facilitated mobile particles and the total number 
of mobile particles at a given time step, provided that they were immobile during the preceding 
time interval At (see text for details). The strain amplitudes from bottom to top are 70 = 
0.02, 0.03, 0.04, 0.05, 0.06. The inset shows the ratio of dynamically facilitated particles and the 
total number of particles that become mobile during the time interval At, given that they were 
immobile during the previous time interval At. The color code for 70 is the same. 
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